import numpy as np
from scipy.optimize import least_squares
import matplotlib.pyplot as plt

# 读取数据文件
with open('智美云科娄底移动光缆在线监测_CH002_2024-12-17_09_11_34-trace.dat', 'r') as file:
    data = file.read()

# 解析数据到 x 和 y 列表中
lines = data.strip().split('\n')
x = []
y = []

for line in lines:
    cols = line.split()
    if len(cols) == 2:  # 确保这一行确实有两个值
        x.append(float(cols[0]))
        y.append(float(cols[1]))

x = np.array(x)
y = np.array(y)

# 使用numpy.polyfit进行多项式拟合，这里我们假设一个二次多项式模型
degree = 6  # 选择适当的多项式阶数
coefficients = np.polyfit(x, y, degree)

# 计算拟合值
y_fit = np.polyval(coefficients, x)

# 计算残差
residuals = y - y_fit

# 定义残差阈值
threshold = 12.0  # 你可以根据实际情况调整这个值

# 找到所有残差绝对值大于阈值的索引
large_residual_indices = np.where(np.abs(residuals) > threshold)[0]

print("Points with residuals larger than the threshold:")
for idx in large_residual_indices:
    print(f"x = {x[idx]}, observed y = {y[idx]}, fitted y = {y_fit[idx]}, Residual = {residuals[idx]}")

# 绘制原始数据、拟合曲线及其残差
plt.figure(figsize=(10, 5))

# 原始数据和拟合曲线
plt.subplot(1, 2, 1)
plt.plot(x, y, 'o', label='Original Data')
plt.plot(x, y_fit, '-', label='Fitted Curve')

# 标记出所有残差大于阈值的点
if len(large_residual_indices) > 0:
    plt.plot(x[large_residual_indices], y[large_residual_indices], 'ro', label='Large Residual Points')
plt.title('Data and Fitted Curve')
plt.legend()

# 残差图
plt.subplot(1, 2, 2)
plt.plot(x, residuals, 'o-', label='Residuals')

# 标记出所有残差大于阈值的点
if len(large_residual_indices) > 0:
    plt.plot(x[large_residual_indices], residuals[large_residual_indices], 'ro', label='Large Residual Points')
plt.axhline(y=threshold, color='g', linestyle='--', label=f'Threshold (+{threshold})')
plt.axhline(y=-threshold, color='g', linestyle='--', label=f'Threshold (-{threshold})')
plt.title('Residuals of the Fit')
plt.legend()

plt.tight_layout()
plt.show()